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Abstract 

We use transition path sampling to study evaporation in the SPC/E model of Uquid water. 
Based on thousands of evaporation trajectories, we characterize the members of the transi- 
tion state ensemble (TSE), which exhibit a liquid-vapor interface with predominantly negative 
mean curvature at the site of evaporation. We also find that after evaporation is complete, the 
distributions of translational and angular momenta of the evaporated water are Maxwellian 
with a temperature equal to that of the Uquid. To characterize the evaporation trajectories in 
their entirety, we find that it suffices to project them onto just two coordinates: the distance 
of the evaporating molecule to the instantaneous liquid-vapor interface, and the velocity of 
the water along the average interface normal. In this projected space, we find that the TSE is 
well-captured by a simple model of ballistic escape from a deep potential well, with no addi- 
tional barrier to evaporation beyond the cohesive strength of the Uquid. Equivalently, they are 
consistent with a near-unity probability for a water molecule impinging upon a liquid droplet 
to condense. These results agree with previous simulations and with some, but not all, recent 
experiments. 

Keywords: Molecular dynamics, rare events, liquid- vapor interface, mean curvature, tran- 
sitition state ensemble, free energy profile 
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Introduction 



In a sample of water at equilibrium with its vapor, the rate of evaporation is equal to the rate 
of condensation. During condensation, not every gas molecule that impinges on a liquid surface 
necessarily sticks. The fraction that does stick is known as the uptake coefficient, 7, and by mi- 
croscopic reversibility, yean also be used to characterize evaporation.^ Any deviation of 7 from 1 
signals some impediment to evaporation (or condensation) beyond the mere cohesive strength of 
the liquid. Measurements of 7 have ranged from about 0.001 to 1 over the past century,!^ but over 
the last decade,*^ they have been converging to the range of 0.1 to 1. Li and coworkers'^ measured 
uptake of isotopically labeled water vapor in a train of water droplets to obtain 7 = 0. 17 ± 0.03 at 
280 K, which increases with decreasing temperature. A similar result, 7= 0.15 ±0.01 at 282.5 K, 
was obtained by Zientara and coworkers'^ from observations of freely evaporating water droplets 
levitated in an electrodynamic trap. Winkler and coworkers,!^ on the other hand, measured droplet 
growth in cloud chambers and claim to exclude values of 7 < 0.4 for temperatures below 290 K. 
Their data is, in fact, consistent with 7=1 for temperatures ranging from 250 K to 290 K. Exper- 
iments done in the Saykally and Cohen groups, ^^I^HIo! which measure the drop in temperature as 
water from a droplet in a droplet train evaporates into vacuum, indicate that 7 = 0.62 ± 0.09 with 
little or no temperature dependence between 245 K and 298 K. 

The experimental uncertainty makes it unclear whether or not there is a small barrier to evap- 
oration. To address this uncertainty, we were motivated to carry out a detailed simulation study 
of evaporation using transition path sampling (TPS), a rare-event sampling technique that can 
produce a statistically representative collection of short evaporation trajectories with Boltzmann- 
distributed (NVT) initial conditions and energy-conserving (NVE) dynamics. Roughly speaking, 
at 300 K, one water molecule evaporates from a 1 nm^ patch of a liquid- vapor interface every 10 ns, 
which motivates using a rare-event sampling technique. Other approaches could also be used to 
study evaporation. For example, in Ref. [I2l a single long simulation of a small water droplet was 
performed at 350 K, resulting in 70 evaporation events. Another complementary approach is to 
study condensation probabilities, since condensation is not rare at all.!^^^ A full discussion of 
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the relationship between evaporation and condensation trajectories is given in the Appendix. The 
chief advantages of our approach are that we do not need to introduce the approximation that the 
velocities and angular momenta of the evaporated water molecule are Boltzmann-distributed, with 
a temperature equal to that of the liquid, and that we are able to generate a large number of evap- 
oration trajectories (about 5000), which we can characterize statistically instead of anecdotally. 
Further, the framework for analyzing TPS simulations can be used to obtain novel insight into 
evaporation kinetics. 

Methods 

Throughout, we run simulations of liquid water with using the SPC/E model of wa- 

ter .^^^ Lennard- Jones interactions are truncated and shifted at a distance of 10 A. Electrostatic in- 
teractions are calculated using the particle -particle particle-mesh (PPPM) method,'^ with a rela- 
tive error of lO^'^. The bond and angle constraints of the water molecule are enforced using the 
SETTLE algorithm's to guarantee that trajectories are time-reversible. A timestep of 2fs is used 
throughout. In simulations where we fix temperature, we use a Langevin thermostat with a time 
constant of 2ps. 

We use the SPC/E model of water because it adequately captures a broad swath of liquid wa- 
ter's properties. With respect to bulk properties at 298 K, its radial distribution function is quite 
accurate,'^ its density is within 1% of experiment,'^ its compressibility'^ of 4.1 x 10^^" Pa^^ is 
close to the experimental value of 4.5 x lO^^^Pa^^ and its dielectric constantl^of 70 compares 
well with the experimental value of 78.2. The properties of its vapor-liquid transition are also 
quite good: the model is explicitly parametrized to reproduce the experimental enthalpy of vapor- 
ization,'i^its liquid-vapor surface tension is within about 10% of the experimental value,'^^^^!] and 
its vapor pressure is within a factor of 2 of the experimental value.*^ With regards to transport prop- 
erties, its self-diffusion coefficientP^l of about 2.8 x lO^^cm^/s compares well the experimentally 
measured value of 2.3 x 10^ cm^/s. These properties lead us to believe that the SPC/E model cap- 
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tures sufficient water-like behavior to be useful in our study. The model is not polarizable, but its 
parametrization accounts implicitly for polarization in the bulk and results in a semiquantitatively 
correct description of the liquid-vapor interface. Of course, it is impossible to obtain arbitrarily 
precise quantitative agreement with experiments using SPC/E or any other classical model of wa- 
ter. However, the consistency of this model with general measures of liquid-vapor coexistence, 
interfacial energetics, and molecular fluctuation amplitudes and time scales gives us confidence in 
its qualitative predictions about the molecular dynamics of water evaporation. 

Transition path sampling^^is used to generate nearly 5000 independent evaporation trajectories 
of length Bps, which is long enough to avoid spurious biases (see Supplementary Information). 
Trajectories are constrained to start in a basin A in phase space, corresponding to a condensed state, 
and end in a basin B, corresponding to an evaporated state. In the analysis below, we consider only 
the trajectories for which the system enters basin B after at least 2ps to avoid any biases towards 
unusually short evaporation trajectories. 

Our explicit definitions of basins A and B are as follows. Basin A consists of all configurations 
where every water that is not hydrogen-bonded to any other water is at most 4 A away from the 
nearest water (the "position of a water" means the position of its oxygen atom, unless otherwise 
stated). Following Ref. |28l two waters are considered hydrogen bonded if the distance between 
their oxygen atoms is below 3.5 A and the angle between the OH bond of the donor and the line 
connecting the two oxygen atoms is below 30°. For our purposes, any other reasonable definition 
of a hydrogen bond should yield nearly identical results. Basin B is defined as all configurations 
of the system where there is exactly one water molecule with no hydrogen-bonding partner that 
is more than 8 A away from its nearest neighbor. The distance cutoff used in defining basin A is 
motivated by the extremely low likelihood for a water in bulk to be that isolated. We comment on 
our choice for the cutoff for basin B below. 

Transition path sampling is essentially a biased random walk in trajectory space. The initial tra- 
jectory of this walk is prepared as follows. We place water molecules in a crystalline arrangement 
in a 30 X 30 X 30 A-^ periodic box so that the density of water molecules matches the bulk density of 
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water (33.3 waters per nm^, equivalent to 0.997 g/ml), i.e., 900 waters in total, and equilibrate this 
system at 300 K for 50ps. Next, we enlarge the box to three times its size along the z-dimension, 
and equilibrate the resulting system for another 50ps. At this stage, we have a 30 x 30 x 90 

D 

periodic box containing a 30 A-thick slab of water parallel to the xy-plane. We then add a water 
molecule about 15 A above the top of this slab with a random, thermal velocity. An example of 
the system at this stage is shown in Figure [T| Next, we evolve the system without a thermostat to 
yield a 3 to 9ps-long trajectory. If in this time, the water molecule does not condense (i.e., enter 
basin A), we discard this initial trajectory and start over. Otherwise, we time-reverse the 3ps stretch 
of the trajectory immediately preceding condensation, and use this reversed evaporation trajectory 
to seed the TPS random walk. We have verified that condensation fails to occur about 50% of 
the time, and that in all cases is due to the water molecule having initial total momentum with a 
positive z-component, so that the molecule moves away from the water slab during the trajectory. 
We have not observed any initial trajectory with a water molecule initially headed towards the wa- 
ter slab and not condensing, a fact that is consistent with a sticking coefficient 7 of nearly 1, as 
observed in previous similar simulations. 1^^^ 

The TPS random walk is performed as follows. At every step, we choose to make a shifting 
move 90% of the time, and a shooting move 10% of the time, reflecting the low cost of shifting 
versus shooting. In a shifting move, we shift the trajectory forwards or backwards by a time At 
uniformly distributed between — 1 and 1 ps. Shooting moves are performed as in the appendices of 
Refs. IS^fand 130, Briefly, the 3A^-dimensional vector of velocities weighted by the square root of 
the atomic masses is rotated slightly, then projected down to a hyperplane to enforce the constraints 
on velocities imposed by the fixed bonds and angles of the water molecules. The kinetic energy 
of the system is then perturbed slightly. Generation and acceptance probabilities for this move are 
chosen to satisfy detailed balance, and the magnitude of the perturbations is chosen to yield an 
approximately 40% acceptance rate. 

For each set of initial conditions, we performed between 10,000 and 20,000 TPS steps, record- 
ing a trajectory every 100 TPS steps. Each recorded trajectory is reasonably independent of the 
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Figure 1 : Snapshot of setup used to study water evaporation. 

previous one, and the first 20 recorded trajectories, which form the equilibration part of the random 
walk, are discarded. To further improve the sampling, we repeated the entire procedure outlined 
in this section about 40 times. The final outcome of this exercise is a set of 4696 mostly un- 
correlated evaporation trajectories, with initial conditions drawn from a canonical ensemble at 
temperature 300K and evolved in time with energy-conserving Newtonian dynamics. 

Our procedure induces a bias for evaporation trajectories where a single water molecule comes 
off the liquid. This bias arises from our definition of basin B for the TPS random walk. Before 
settling on this definition, we explored the possibility of events where dimers or larger aggregates 
of water evaporate as a unit, by using a more generous but cumbersome definition of basin B. 
Specifically, a configuration was in basin B if it contained exactly two separate clusters of waters, 
in each of which every water was close to some other water in the cluster. By observing the 
evaporation events in these preliminary simulations, we convinced ourselves that out of the rare 
events in which evaporation occurs, those involving more than one water were far rarer still, so we 
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neglected this possibility in our final simulations in favor of using a simpler definition of basin B. 

In analyzing the evaporation trajectories, it is useful to locate the hquid-vapor interface at any 
instant in time, for which we use the method of Ref . '3T, Briefly, we map a given configuration of 
water oxygen atoms {r,} onto a smooth density field p(r) defined by the relation 

p(r) = £<^)(|r-r,-|), (1) 

where N is the number of water molecules, and ^(r) is a Gaussian-hke smoothing function of 
width ^ = 2.5 A (see Appendix). The instantaneous liquid-vapor interface is then defined implicitly 
as the set of points {s} that satisfy 

p(s) = (l/2)p^, (2) 

where pi is the bulk density of liquid water. 

After locating the liquid- vapor interface, we follow Ref. 31 in defining a perpendicular dis- 
tance a from a probe water molecule at r to the interface, as illustrated in Figure [2j First, we 
locate the point s on the interface closest to r, and calculate the vapor-pointing normal vector to 
the interface there, n. Then a is defined as the distance from r to s projected along the n direction, 

a = n-(r — s). (3) 

In Ref. |3T1 this distance was denoted by a* . 

An ambiguity arises about whether the probe molecule at r should or should not be included 
when calculating the position of the liquid-vapor interface. Generally, we exclude it when calcu- 
lating a. To discuss the consequences of this choice, we define a' analogously to a, but with the 
probe water molecule included in the definition of the liquid-vapor interface. 
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(a) (b) 

Figure 2: Definition of water-to-surface distance a (or a'). The plane shown is the tangent plane to 
the liquid- vapor interface at the point closest to the probe water molecule. In (a), the probe water 
is excluded from the definition of the interface; in (b), it is included. These snapshots illustrate a 
typical interfacial deformation that accompanies an evaporation event. The system is depicted at 
its transition state. 

Results 

Evaporation correlates with negative mean curvature 

We first focus on the molecular details of the transition states of the evaporation trajectories. Or- 
dinarily transition states are identified using committor functions. The committor, pb(x), of a 
spatial configuration x is defined as the fraction of short trajectories that start at x with random 
thermal velocities, and finish in basin B. At most points in a transition path, this function is either 
or 1, with a quick crossover around the configurations that dominate the dynamical bottleneck 
between A and B. Thus, a pragmatic definition of a transition state along a trajectory is the point 
where j9g(x) = 0.5. 

Implicit in the above definition of the committor function is the assumption that momenta are 
not important in characterizing transition states. In a dense system, this assumption is generally 
true, since the velocity of any particle decorrelates rapidly, usually within 1 ps.'^ When examin- 
ing evaporation, the assumption breaks down, since the velocity of an evaporated water molecule 
decorrelates over much longer timescales. The clearest manifestation of the problem is that the 
standard definition of PB{t) leads to psit) ~ 0.5 for a configuration containing a single, clearly 
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evaporated water molecule, since the water can likely recondense if its rethermalized velocity 
points towards the liquid slab. 

As a compromise, we have chosen to redefine the committor function to include the z-component 
of the velocity of the evaporated water molecule. Strictly speaking, it's impossible to tell which 
water molecule is "the evaporated molecule" in an arbitrary configuration, but this is not a problem 
for identifying transition states along a transition path. Figure [3] illustrates the typical behavior of 
PB{t) = pBi'^it) ,v^z^^ {t)) defined in this way, estimated by spawning 10 short trajectories at every 
time point. 




Figure 3: Estimated committor, pgit), sampled at 20 fs intervals for several evaporation trajecto- 
ries. The error in pg in the region around 0.5 is about 0.15. To obtain this estimate, all velocities 
but the z-component of the evaporated water's velocity are randomized independently 10 times, 
after which a short trajectory is evolved forwards in time for up to 5ps until the systems enters 
either basin A or basin B. For clarity, individual committors are slightly displaced vertically. 

We have defined the transition states as the configuration at a time tc equal to the mean of 
the first time for which pB{t) exceeds 0.4 and the first time for which it exceeds 0.6. The exact 
value of tc is not very sensitive to the chosen cutoffs, as long as they are reasonable. The set of all 
configurations of the evaporation trajectories at their respective times tc comprises the transition 
state ensemble (TSE). 
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In many condensed-phase phenomena, collective coordinates are key. Positions of individual 
atoms in the TSE do not by themselves appear particularly remarkable or extraordinary. Visual 
inspection confirms this state of affairs in this particular case. Instead, it is essential to characterize 
the members of the TSE by looking for statistical trends in a few collective coordinates. Here, we 
focus on the instantaneous liquid-vapor surface. Let s be the point on this surface that is closest 
to the evaporating water molecule at any given time. The mean curvature, H, of the surface at s 
serves as a concise characterization of collective fluctuations of water molecules at the liquid-vapor 
surface. The mean curvature is defined 



H=''A^. (4) 



where ki and ^2 are the principal curvatures at s. The magnitude of a principal curvature is the 
reciprocals of the principal radius of curvature, and its sign specifies whether the surface curves 
towards (positive) or away (negative) from the normal direction along the corresponding principal 
direction. The mean curvature characterizes the change in surface area upon infinitesimal defor- 
mation of the surface, so it can be interpreted as a local characterization of the force of surface 
tension on the liquid-vapor surface. In particular, a deformation along the normal direction by an 
infinitesimal distance £ changes the area element dA as^l 



dA^ (l-2e//)dA. (5) 

To establish a baseline, we first calculate the distribution of // as a function of the height a of 
a probe water molecule from the liquid- vapor interface. Figure |4] shows the results as a joint free 
energy for H and a (respectively, H' and a' if the probe water molecule is included in the definition 
of the liquid-vapor interface), calculated using umbrella sampling as described in the Appendix. 
At very low and very high a, only a trivial bias in H is seen as a function of a, resulting from the 
nearest point on the surface being preferentially one where the surface is bending most towards the 
probe water molecule. However, an evident additional bias towards negative mean curvature can 
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be seen for a just above the surface, indicating that a water molecule suspended there significantly 
deforms the surface below it. Figure [2] shows an example of this kind of deformation in one of the 
harvested evaporation trajectories. 
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Figure 4: Free energy for height a of a probe water molecule and the mean curvature, H, at the 
nearest point on the liquid-vapor surface (respectively a' and H' if the probe water molecule is 
included in the definition of this surface). Contours are spaced at Ik^T. 

Figure |5] overlays the transition states of the evaporation trajectories on the free energies of Fig- 
ure |4j To a certain extent, the transition states exhibit some of the bias towards negative curvature 
that can be seen in the equilibrium free energies. The bias is slight when the probe molecule is 
not included in the definition of the liquid-vapor surface, but is clearer when the probe molecule 
is included. The definition of a liquid-vapor interface during the evaporation process is somewhat 
ambiguous, and we regard full inclusion and full exclusion as the two limiting extremes for a 
suitable definition. Since the bias towards negative curvature is present in both cases, our finding 
should be robust with respect to reasonable changes in the definition of the interface. 

As discussed below, the preponderance of negative -mean-curvature liquid-vapor interfaces 
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Figure 5: Representative transition states of evaporation trajectories (red) projected onto the 
H and a coordinates. The free energies of these coordinates are shown for comparison. Labels 
as in Figure |4j 
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does not correspond to an entropic barrier to evaporation, but instead is a molecular manifesta- 
tion of the cohesive strength of the liquid. Nevertheless, we anticipate that external influences 
might be used to alter the microscopic details we describe, and so may perhaps be used to exert 
control over evaporation. Additionally, our characterization establishes a baseline for understand- 
ing evaporation under different conditions where barriers are observed in simulations, such as at 
higher temperatures'^ or in the presence of surfactants.'^ 

Post-evaporation momenta are Boltzmann-distributed 

We now examine the center-of-mass velocities and angular momenta at the end of each trajectory. 
In all of the following results, we first estimate the value of each observable independently in each 
TPS run, and then report the mean of these values, with an error bar estimated as the standard error 
of the mean. 

Figure[6]^a) shows the distributions of the component of the evaporated water molecule's center- 
of-mass velocity along a direction perpendicular to z, measured at the end of an evaporation trajec- 
tory. Figure |6jb) shows the analogous distribution of the components of angular momenta along the 
principal axes of inertia of the evaporating water molecule. Both sets of distributions are consistent 
with Boltzmann statistics at temperature T = 300 K. 

The component of the velocity along the z direction, v^, has a more interesting distribution, 
shown in Figure ^c). We enforce the constraint that water molecules first enter basin B with a 
positive by flipping trajectories where this is not the case. Hence, no water molecules should 
have negative at the end of the evaporation trajectory if the definition of basin B were sufficiently 
strict. In practice, the definition of basin B used here does not perfectly discriminate between the 
evaporated states and states where recondensation will occur. Since the trajectories examined here 
are finite, a trajectory where the system that transiently enters B before recondensing may appear as 
an evaporation event, but with Vj; < at the end of the trajectory. Only about 1% of our trajectories 
exhibit this problem, which can in principle be mitigated by using longer trajectories and a stricter 
definition of basin B. 
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Figure 6: Distribution of (a) the component of center-of-mass velocity perpendicular to z; (b) the 
components of angular momentum of the evaporated water along the principal axes of inertia, 
measured at the end of an evaporation trajectory (symbols); and (c) the z-component of center-of- 
mass velocity. For (a) and (b), the relevant Boltzmann distributions at temperature T = 300 K are 
also shown (dashed lines). For (c), the expected result for thermal ideal gas particles evaporating 
from a deep, barrierless potential well is shown (Equation ([10]), dashed line). 
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The expected distribution of for positive can be deduced from a simple model (Figure |7]) 
of thermal ideal gas particles evaporating from a deep, barrierless potential well of depth Ai7. 
Particles inside the well have a thermal distribution of velocities, P{vi), given by 

/'(v,)ocexp(^-^/3mv2^. (6) 

A particle with initial velocity velocity v, can only escape the well if is above a threshold velocity, 
Vt, given by jtnvj = AU. Were there a barrier, this threshold velocity would be higher, but the 
remainder of this discussion would carry through unchanged. The final velocity of this particle, Vf, 
is determined by conservation of energy, independent of the details of any intermediate barrier: 

1 9 1 9 

-mvt = -mv} + AU. (7) 

This equation relates the distributions of initial and final velocities, P{vi) and P{vf) respectively, 
after correcting for the fact that for finite trajectories, high initial velocities are overrepresented 
by a factor of |v,|, as there are proportionally more possible starting positions compatible with the 
particle being outside the well at the end of the trajectory. The exact relationship is 

P(v/)dv/ocP(v,-)|v,|dv„ (8) 



so 



P(vf) ^ P(v,)|v,|^ . exp (-ipmv?) -j^. (9) 



Since the denominator in the last fraction is equal to |v/|, we have 

,^,^^^^•Sr^■M-^-^^"'■})' "/><'• (10) 
0, Vf < 0. 

Were there a barrier of height B to evaporation, the threshold v / above would be ^JlB/m instead 
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of 0, but the functional form would remain unchanged. 



U{z) 
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Figure 7: Ideal gas particles at the bottom of a deep, barrierless potential well have a Boltzmann 
distribution of velocities. Only a fraction of particles have enough energy to escape the well. 
After evaporating, but before thermalizing outside the well, the distribution of velocities of these 



particles is given by Equation (10) 



While Equation ([TO]) was derived for an ideal gas of thermal particles escaping from a deep, bar- 
rierless potential well, it also follows more generally from considerations of time reversibility (see 
Appendix) and it describes the observed distribution of Vy for evaporating molecules surprisingly 
well (dashed line in Figure [6];c)). In particular, low- velocity particles are not underrepresented, 
which is consistent with there being no barrier to evaporation. A similar velocity distribution has 
been reported in simulations of argon evaporation, which can be observed straightforwardly with- 
out special sampling techniques like TPS.I^ 



Potential of mean force for removing a water molecule from bulk is barrier- 
less 

Figure [S] shows the free energy, F{a), of an arbitrary water molecule in our system as a function of 
the perpendicular distance to the instantaneous hquid-vapor interface a, calculated using umbrella 
sampling (see Appendix). Such a free energy profile is a reversible work or a potential of mean 
force surface (i.e., its negative gradient is equal to the mean force experienced by a water molecule 
along the a coordinated^). The essential feature of this free energy is that it is barrierless. Apart 
from density layering in the bulk,^^ manifested as oscillations in F{a) for a < OA, the bulk liquid 
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simply sets up a deep potential well for any individual water molecule, and a molecule in the vapor 
can simply roll downhill into this well. While the absence of a barrier along the a coordinate 
does not preclude the existence of barriers along other coordinates, we demonstrate below that the 
transition states of the evaporation trajectories are consistent with a describing the majority of the 
evaporation reaction coordinate. 
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Figure 8: Free energy for a single water molecule at a perpendicular height a from the liquid-vapor 
interface defined by the remaining water molecules. The red lines are the free energies of the stable 
liquid and vapor phases, and are guides to the eye. The biasing potentials used extend to a = 7 A, 
so the apparent downturn at a = 8 A is not statistically significant. 

The depth of the well in F{a), denoted by AF*, quantifies the cohesiveness of the liquid with 
respect to the vapor. Indeed, if we regard a single water molecule as an independent particle 
moving in the potential well F{a), then the relative density of this particle in the liquid, pf, with 
respect to that in the vapor, pg, is given by 

Pg = P(e-^^''\ (11) 

We estimate from Figure [8] a value of AF* of 11.5 ±0.2^3 7"- This compares favorably with the 
value of ll.^kBT obtained by setting pg = Py^cg/k^T and using the computed value of /Vap for 
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SPC/E water at a temperature of 300 K and pressure of latm-f^^For real water, the analogous 
calculation yields AF* = 10.5 ksT. 

The range of F{a) also characterizes the effective range of attraction between a molecule in 
the vapor and the bulk slab, just under 8 A. It is this range that motivates the definition of basin B 
described in the Methods section. Different models of water will have slightly different ranges of 
attraction, but we do not expect discrepancies in the qualitative behavior of F{a). 

Others have calculated a similar potential of mean force, but with respect to the distance from 
the Gibbs dividing surface instead of the instantaneous liquid surface, so that the details of the 
potential are masked by the capillary wave fluctuations of the liquid-vapor interface. Nevertheless, 
their results for the SPC/E water model and for a polarizable water model due to Dang and 
Chang — are broadly similar to each other and to our own results. 

Transition states are consistent with diffusion out of a deep well 




Figure 9: Evaporation trajectory traces projected onto variables vf^^'^'' and a (black lines). The 
transition state of each trajectory, identified as described in the text, is highlighted by a green dot. 
Red line: the expected transition state ensemble for a coarse model of ballistic escape from a 
potential shaped as in Figure [81 given by Equation ([T2[). 
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Figure |9] depicts traces of many evaporation trajectories projected onto the two coordinates 
^(evap) ^.j^j^ transition state of each trajectory highlighted in green. Unlike similar traces 

onto many other pairs of coordinates (not shown), there is a definite correlation between the dis- 
tance of the evaporated water from the liquid-vapor interface and its speed in the z direction. We 
can partially rationalize this dependence by conceiving of the free energy along a (Figure [8]) as 
an actual potential energy well, and approximating the velocity along the a direction with v^^^''^^ 
If evaporation were a ballistic escape from this well, then the transition states would satisfy the 
condition 

^m(v('=^^P))2 = F(a). (12) 

The points satisfying this relation are shown as a thick red line in Figure |9j Despite the evident 
crudeness of the model, the transition states clearly cluster around the line of Eq. ( [12] ). 



Discussion 

We have examined the process of evaporation of SPC/E water in detail, and all the evidence sug- 
gests that there is no barrier to evaporation in this model. In other words, to evaporate, a water 
molecule near the surface only needs to spontaneously acquire enough kinetic energy in the di- 
rection of the liquid-vapor interface normal. This view is consistent with the distribution of 
for the final velocities (Figure [6];c)), the fact that the potential of mean force along a coordinate a 
perpendicular to the liquid-vapor surface is barrierless (Figure [8]) and the fact that the transition 
states cluster around values of and a that have a threshold amount of energy to escape from the 
potential well set up by the remainder of the bulk (Figure |9]). It is difficult to imagine evaporation 
to be a mildly activated process and still be consistent with these three pieces of evidence. 

Our results are consistent with the near-unit condensation coefficient measured in simulations 
in Refs. [121 [14] and [121 but is in apparent contradiction with the most recent experimental re- 
sults,'^'^^^ which suggest a barrier of around —^87^111(7) ~ 0.5 ^bT"- The other experimental re- 
sults cited in the introduction suggest anything from the absence of a barrier to a barrier of up 
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to l.9k^T. Excluding the possibility that water molecules evaporate as dimers, which would im- 
ply that an appreciable fraction of water molecules in the vapor as dimerized (and recall that our 
preliminary transition path sampling showed that there is not a significant fraction of SPC/E water 
molecules that evaporate or condense as dimers or as larger clusters), such large barriers should be 
clearly evident in direct simulations of water condensation, but they are conspicuously absent.'^^'^ 

The general lack of consensus between experiments'^ makes it unclear whether or not our 
result of apparent unit evaporation coefficient agrees with reality, or if it is an artifact of our sim- 
ulations. In particular, it could be that there is indeed a barrier to evaporation and we cannot 
capture it, if that barrier were due to fundamentally quantum effects. By construction, these ef- 
fects are beyond the scope of the classical molecular dynamics simulations used here. Important 
quantum effects are plausible because librational motions of water are strongly quantized: their 
typical wavenumbers, around 500 cm^, are comparable to the thermal wavenumber atT = 300 K, 
around 200 cm ^ More sophisticated simulation techniques can incorporate many quantum effects 
at reasonable cost. A notable exception would be dynamical quantum coherence, 1^ for which a sig- 
nificant role would be surprising for intermolecular motions in a strongly dissipitating system like 
liquid water. If quantum effects were limited to quantum dispersion and simple tunneling behav- 
ior, for instance, one could explore the consequences of quantum uncertainty using ring-polymer 
molecular dynamics.'^ However, our firm expectation is that these more sophisticated simulations 
will produce results that agree with those presented here, since generally, tunneling and dispersion 
tend to lower effective barriers with respect to classical expectations, not increase them. Moreover, 
any account of such quantum effects playing a dominant role would have to be compatible with 
the observation^ that the evaporation coefficient of D2O is equal, within errors, to that of H2O. 

Another possible source of discrepancy is our use of the SPC/E model of water, and in particu- 
lar, its lack of polarizability, which might result in a qualitatively inaccurate description of events 
at the liquid-vapor interface. However, the agreement of its surface tension to the experimental 
value (within about 10%) suggests that the SPC/E model's parameters implicitly capture enough 
detail about polarization to describe the general mechanistic behavior of the liquid- vapor interface. 
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Moreover, the addition of polarizability would likely reduce, not enhance, any barriers to evapo- 
ration and/or condensation, since polarization induces an additional attractive force between the 
liquid and a vapor molecule that is relatively long-ranged. Significantly, a previous study of direct 
condensation that used the P0L3 model of water, which is polarizable, is consistent with 7 ~ 1 , 
i.e., barrierless evaporation. 

Finally, extracting the evaporation coefficient from experiments involves some interpretation 
and extrapolation, so it is conceivable that the quoted results may be skewed by systematic errors 
that have not been accounted for. For example, Morita et al.'^ have previously argued that Li et 
al's low reported evaporation coefficient!^ may actually be compatible with a value in the range of 
0.2 to 1 once the effects of fluid flow on the gas surrounding their water droplet train are corrected 
for. As for the more recent experiments of Refs. [3l [H |9]and [TOl these rely on a linear extrapolation 
of van't Hoff behavior of the Raman spectrum of water down to supercooled temperatures in order 
to measure the temperature of evaporating water droplets. Recent Raman spectra of magnetically 
trapped supercooled droplets, however, show that this extrapolation may not be accurate. This 
suggests that the observed deviation from unit evaporation coefficient may also be in part due to 
shortcomings in the calibration step of the experiments. A systematic error of 2 % in absolute 
temperature in the experiments (equal to about 10 % in the temperature change during the course 
of the measurements) would be sufficient to account for the discrepancy between the experiments 
and our calculations. 
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Appendix 

Time reversibility and evaporation vs. condensation 

The observables measured using TPS for evaporation can be related to those measured in simu- 
lations of condensation. Let hA{x) and hB{x) be indicator functions of basins A and B. They are 
equal to 1 if the phase space point x is in the respective basin, and otherwise. The Boltzmann dis- 
tribution, which specifies the initial conditions for our evaporation trajectories, is denoted by p{x). 
The quantity Pj [x — y] is the probability density that a trajectory of length T has its endpoint in 
the vicinity of y, given that it started at x. For the energy-conserving dynamics that we use in the 
text, 

PT[x-^y] = 5\y-UT{x)i (13) 

where Ut{x) is the time evolution operator over a time T. 

The expectation of an observable W measured at the endpoint b of an evaporation trajectory of 
the kind sampled by TPS is given by'^ 

_ Jdb JdahA{a)p{a)PT[a ^ 
^^I'^^/evap JdbJdahAia)p{a)PT[a^b]hB{b) ' ^ ^ 

Conversely, the expectation of ^ measured at the beginning point a of a condensation trajectory 



23 



can be defined as follows: 



From these definitions, the following relation follows immediately: 

{JdahA{a)PT[a^b])^^^^ 

Any configuration x can be mapped onto its time-reversed counterpart, which we denote x, by 
inverting the direction of all the particle momenta. For time-reversible dynamics, such as that used 
in the text, we have 

Pria ^ b] = Prib ^ d], (17) 
and further, for energy-conserving dynamics, if a and b are in the same trajectory, then 

p{a)=p{b). (18) 



With these relation, we can rewrite Eq. ( 16 ) in a more usable form. 



^ {^{b)JdaPT[b^a]hAia))^^^^ 
{JdaPT[b-^a]hA{a))^^^^ 



(19) 
(20) 



In the second equation, we have renamed the integration variables a and b, and exploited that 
= hA{a) and hsib) = hgib). 

Equation ( [20| ) tells us that averages over TPS trajectories are equivalent to time-reversed av- 
erages over trajectories that start in B and end in basin A after time T . A priori, there is no 
requirement that the water that is condensing have an initial velocity that is directed towards the 
liquid slab, though trajectories that do not satisfy this condition are very unlikely to end in basin A. 

A subtle point about Equation pO] ) is that the conditional factor / daPr — )■ a\hA{a) cannot be 
approximated as 1/2 for large T. Indeed, basin B is potentially unbounded, so no matter how large 
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a r is chosen, there will be configurations in B with an initial velocity of the isolated water is too 
low for the system to escape basin B in time T. There are two potential solutions to this problem. 
One is to make basin B finite. Alternately, and more revealingly, one can model the consequence 
of the unboundedness of basin B, as we do below. 

For concreteness, we consider a simpler definition of B than the one used in the text, which 
is easier to analyze and allows us to make the connection between the discussion here and kinetic 
rate theory.!^ Let z{b) be the z-coordinate of the evaporated water molecule's center of mass, and 
let Vz{b) be the corresponding component of the velocity. The simpler basin B consists of all con- 
figurations b for which z{b) > z*. With this definition, we can make the following approximation: 

Prib ^ a]hA{a) ^ e[\v,ib)\T - {zib) - z^^ib* ^ a], (21) 

with T ^ r a small, fixed time and b* the point along the trajectory starting at b where z is first equal 
to z*. In other words, the probability for a configuration b to end in basin A is mostly determined 
by whether T is long enough to get to the boundary of B, and then a kinetic factor that's virtually 
independent of T. We also assume that W{b) is independent of z{b), so we can replace ^(Z?) by 
) . Since the mapping from b to b* is area preserving, we have 



^ {nb)\vz{b)mb)-z*] IdaP,[b ^ a]hA{a)),,,, 
{\v,{bmz{b)-z*]!daP,[b^a]hM),,,, ' 

In comparison, the transmission coefficient for a reaction from 5 to A after a transient time T is 
given byl^ 

^ . . _ {\vmS[z{b)-z*]JdaP,[b^a]hA{a)),,,^ 

(|v.W|5[z(Z.)-z1-(l/2)U, • ^'^^ 

As is normal in reaction rate calculations, this transmission coefficient is almost independent of 
T for values of T greater than molecular timescales but smaller than implied by typical reaction 
rates. Here, those conditions require that 1 ps ^ T ^ 1 ns. In this plateau regime, the transmission 
coefficient is equal to the uptake coefficient, 7. If this coefficient is 1 and z{b) is high enough that 
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an initially evaporating water molecule does not recondense, then we have 



JdaPr[b a\hAia) ^ &[-Vz{b)], (24) 



so that 



{^ ib)\v,ibmz{b)-z*]&[-vm),o ,i 
\vm5m-z*M-v,{b)])coM 



The quantity on the right-hand sides of Equations ([22]) and ( [25] ) is what is directly measured in 
condensation simulations. Obtaining them required several assumptions, all of which are reason- 
able in the context of this paper. However, our treatment here highlights the assumptions explicitly, 
and will be useful in contexts where these assumptions may not apply. 

One simple application of Equations ([22]) and ([25]) is to calculate the distribution of Vy for the 
evaporated water molecules. Substituting ^(Z?) = d[vz{b) —Vf] immediately yields Equation (]9]). 



Choice of density smoothing function (j) (r) 

In the main text, the liquid-vapor interface is defined as an isosurface of the smoothed density 
field p(r), constructed by convoluting the instantaneous water density (a sum of Dirac delta func- 
tions) with a smoothing kernel, ^(r). In Ref. [3TJ the choice for ^(r) was a Gaussian of width ^, 
truncated and shifted at r = 3^. Since our study focuses on the curvature of the liquid-vapor inter- 
face, the discontinuity in first and second derivatives of (r) at the cutoff point is inconvenient. In- 
stead, to ensure that p(r) is sufficiently smooth, we use a 0(r) that results from stitching two cubic 
functions of r at the point r = c, subject to the following conditions: (a) ^(r), ^'{r) and ^"{r) are 
continuous at r = c, {h)(j){3^) = 0, (c) (j)'{0) = (j)'{3^)=0, (d) (j)"{3^)=0, (e) Jq dr4nr^(l){r) = 1. 
These eight conditions uniquely specify ^(r). The stitching point is chosen empirically to be 
c = 2.1^ so that (^(r) closely resembles a Gaussian with standard deviation ^. In Ref. ]3T1 a value 
of = 2.4 A was chosen, which leads to about 7 % of our trajectories having an ambiguous liquid- 
vapor interface at some timestep (i.e., Eq. (|2]) defining more than two liquid-vapor interfaces). We 
have found it convenient to use a slightly higher value, ^ = 2.5 A, whereby the fraction of trajec- 
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tories with ambiguous liquid- vapor interfaces at any timestep drops to about 3%. For simplicity, 
all of these trajectories are discarded in their entirety in the analyses above. 

Umbrella sampling with respect to the position of the liquid-vapor interface 

We have used umbrella sampling to collect statistics on rare configurations of our system where a 
probe water molecule is at a fixed perpendicular distance a (or a') from the instantaneous liquid- 
vapor interface. To do this, we have used the indirect umbrella sampling method (INDUS) that we 
have previously used in different contexts .'^Briefly, we umbrella sample along a different coor- 
dinate that tracks a, use to properly reweight all our samples, then compute histograms 
for a and possibly other variables from these weighted samples. The coordinate we use is the dis- 
tance a from the probe water molecule to the instantaneous liquid-vapor interface directly below 
it. Let h{x,y;r^) be the z-coordinate of the liquid-vapor interface with the given values of x and y, 
which in turn depends on the positions of the N water oxygen atoms. The umbrella potential we 
use is 

Vir^) = |[z„-/^(x„,j„;0 -a]\ (26) 

Here, n is the index of the probe water molecule, with coordinates (xmymZn)- The value of /z(.)c„,j„;r^) 
is defined implicitly by the equation 

p{xn,ynM^„,y„;r^y,r^) = (l/2)p^ (27) 

We henceforth suppress the dependence of it on r^. In a slab of water, there are usually two disjoint 
interfaces at the slab's top and bottom, so this equation has two solutions. For concreteness, we 
always refer to the top interface of the slab. 

To calculate h{xn-,yn) quickly at every timestep, as well as its gradient with respect to particle 
positions, we note that the value of h{xmyn) at one timestep is similar to its value at the next 
timestep. We have thus implemented a parallel Newton-Raphson solver to calculate h{xn,yn), with 
the starting guess at one timestep equal to the value of h{xn,yn) at the previous timestep. In a 
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T O 

typical simulation, convergence to 10^ A occurs after just one or two iterations. 

To calculate the forces implied by the umbrella potential, we need to calculate the gradient of 



Equation ( [26| ) with respect to particle positions. We present explicit expression below, where h and 
its derivatives are evaluated at (jc„ , y,, ) , while p and its derivatives are evaluated at {xn ,yn-,h{xn-,yn))- 
To simplify the calculation, we assume that the tagged particle n is not a water oxygen, and then 



relax this restriction. By taking the total derivative of Equation ( [27] ) with respect to the position of 
oxygen atom z, we find that 

d(p-p^/2) dp dh dp 
dvi dz dr, dr,- 

Hence, 



dh dp /dp 
dr,- dri / dz ' 

The derivative with respect to the position of particle n is obtained similarly, so 



(29) 



dip-p^^dp^dp^^^^ (30a) 

dxn dx dz dXn 

d(p-pe/2) dp dp dh 
VP Ht/ i_ P (30b) 



dyn dy dz dyn 

d{p-pf/2) 

dZn 



0. (30c) 



Hence, 



dh 


dp 


/dp 


djin 


dx/ 


dz' 


dh 


dp 


/dp 


dj„ 


dy/ 


dz' 


dh 

dzn 


= 0. 





(31a) 
(31b) 
(31c) 

If the probe water molecule n is itself included in the definition of the liquid- vapor interface, then 



d/i/dr„ is the sum of the right-hand sides of Equations (|29]) and pT]) 
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Supplementary Information: Length of evaporation trajectories 

In this section, we show that the 3ps length of our TPS trajectories is long enough. 

For each trajectory, let Ia be the latest time for which the system is in basin A, and let ts be 
the latest time for which the system is not in basin B. These times roughly characterize the points 



along the trajectory at which the evaporation event begins and concludes. Figure 10 shows the 
distribution of the time difference tg — tA- Most evaporation events take under Ips, and very few 
take just under 3 ps. Hence, the 3ps trajectory length we chose to use for our TPS sampling is long 
enough. Correcting the distribution of times tg — tA for the bias towards short evaporation events 
owing to their larger number of possible starting times does not change this conclusion. This is 
demonstrated in Figure [TTJ which shows the distributions of times fg. If the TPS trajectory length 
is sufficiently long, then this distribution should rise from zero at small ?g and plateau to a constant 
for tB much larger than the typical time for an A-to-B transition to occur. This is indeed observed. 
Were the TPS trajectory length too short, there would be no plateau region. 
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